library(tidyverse)
## ── Attaching packages ───────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.1.0       ✔ purrr   0.3.2  
## ✔ tibble  2.1.1       ✔ dplyr   0.8.0.1
## ✔ tidyr   0.8.3       ✔ stringr 1.4.0  
## ✔ readr   1.3.1       ✔ forcats 0.4.0
## ── Conflicts ──────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(ggmap)
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
library(ggExtra)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

3 - Create Illustrations

Prepare color palette

library(viridisLite)
colors <- rev(magma(16))
colors_log <- colors[-c(7, 8, 10, 11, 13, 15)]

Load maps

roadmap <- readRDS("../preprocessed_data/roadmap.rds")
toner_lite <- readRDS("../preprocessed_data/toner_lite.rds")
toner_lines <- readRDS("../preprocessed_data/toner_lines.rds")
napa_roadmap <- readRDS("../preprocessed_data/napa_roadmap.rds")
napa_toner_lite <- readRDS("../preprocessed_data/napa_toner_lite.rds")
napa_toner_lines <- readRDS("../preprocessed_data/napa_toner_lines.rds")

ggmap(roadmap)

ggmap(toner_lite)

ggmap(toner_lines)

ggmap(napa_roadmap)

ggmap(napa_toner_lite)

ggmap(napa_toner_lines)

Prepare blank plot

blankPlot <- ggplot()+geom_blank()+ theme_void()

Load directory

directory <- readRDS("../preprocessed_data/directory_loc.rds")

Remove observations that only say “California” as address.

head(sort(table(directory$address), decreasing=TRUE))
## 
##                                    California 
##                                           487 
##        389 4th St E\nSonoma, California 95476 
##                                            55 
##                        Napa, California 94558 
##                                            53 
##                              Napa, California 
##                                            47 
##                  Healdsburg, California 95448 
##                                            44 
## 389 4th Street East\nSonoma, California 95476 
##                                            37
nrow(directory)
## [1] 20472
directory <- directory[directory$address != "California", ]
nrow(directory)
## [1] 19985

Test plot directory - 2015

directory_map_point <- ggmap(roadmap, base_layer = ggplot(subset(directory, year==2015), aes(lon, lat))) +
  coord_cartesian() +
  geom_point(alpha=0.1, stroke=0, position=position_jitter(0.02, 0.02)) +
  geom_density2d()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
directory_map_point
## Warning: Removed 116 rows containing non-finite values (stat_density2d).
## Warning: Removed 116 rows containing missing values (geom_point).

napa_map_point <- ggmap(napa_roadmap, base_layer = ggplot(subset(directory, year==2015), aes(lon, lat))) +
  coord_cartesian() +
  geom_point(alpha=0.1, stroke=0, position=position_jitter(0.01, 0.01)) +
  geom_density2d()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
napa_map_point
## Warning: Removed 1799 rows containing non-finite values (stat_density2d).
## Warning: Removed 1800 rows containing missing values (geom_point).

Load USDA

usda <- readRDS("../preprocessed_data/usda_long.rds")

Test plot USDA - 2015

directory_map_point <- ggmap(roadmap, base_layer = ggplot(subset(usda, year==2015), aes(lon, lat))) +
  coord_cartesian() +
  geom_point(alpha=0.3, stroke=0, position=position_jitter(0.1, 0.1)) +
  geom_density_2d()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
directory_map_point
## Warning: Removed 2 rows containing non-finite values (stat_density2d).
## Warning: Removed 2 rows containing missing values (geom_point).

napa_map_point <- ggmap(napa_roadmap, base_layer = ggplot(subset(usda, year==2015), aes(lon, lat))) +
  coord_cartesian() +
  geom_point(alpha=0.5, position=position_jitter(0.3, 0.3)) +
  geom_density2d()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
napa_map_point
## Warning: Removed 128 rows containing non-finite values (stat_density2d).
## Warning: Removed 138 rows containing missing values (geom_point).

Final plots - directory progression

directory_2015 <- ggmap(roadmap, base_layer = ggplot(subset(directory, year==2015), aes(lon, lat))) +
  coord_cartesian() +
  geom_point(alpha=0.1, stroke=0, position=position_jitter(0.02, 0.02)) +
  geom_density2d()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
xdensity_2015 <- ggplot(subset(directory, year==2015), aes(x=lon)) +
  geom_histogram() +
  scale_x_continuous(limits=c(-126, -113)) +
  labs(y=element_blank(), x=element_blank())

ydensity_2015 <- ggplot(subset(directory, year==2015), aes(x=lat)) +
  geom_histogram() + 
  coord_flip() +
  scale_x_continuous(limits=c(31.5, 41.5)) +
  labs(y=element_blank(), x=element_blank())

distribution_2015 <- gridExtra::arrangeGrob(xdensity_2015, blankPlot, directory_2015, ydensity_2015, 
                                            ncol=2, nrow=2, widths=c(4, 1), heights=c(1, 4))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 115 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 116 rows containing non-finite values (stat_density2d).
## Warning: Removed 116 rows containing missing values (geom_point).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 113 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
grid::grid.draw(distribution_2015)

ggsave("../images/distribution_2015.png", distribution_2015, width=7, height=7)
directory_2017 <- ggmap(roadmap, base_layer = ggplot(subset(directory, year==2017), aes(lon, lat))) +
  coord_cartesian() +
  geom_point(alpha=0.1, stroke=0, position=position_jitter(0.02, 0.02)) +
  geom_density2d()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
xdensity_2017 <- ggplot(subset(directory, year==2017), aes(x=lon)) +
  geom_histogram() +
  scale_x_continuous(limits=c(-126, -113)) +
  labs(y=element_blank(), x=element_blank())

ydensity_2017 <- ggplot(subset(directory, year==2017), aes(x=lat)) +
  geom_histogram() + 
  coord_flip() +
  scale_x_continuous(limits=c(31.5, 41.5)) +
  labs(y=element_blank(), x=element_blank())

distribution_2017 <- gridExtra::arrangeGrob(xdensity_2017, blankPlot, directory_2017, ydensity_2017, 
                                            ncol=2, nrow=2, widths=c(4, 1), heights=c(1, 4))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 121 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 122 rows containing non-finite values (stat_density2d).
## Warning: Removed 122 rows containing missing values (geom_point).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 118 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
grid::grid.draw(distribution_2017)

ggsave("../images/distribution_2017.png", distribution_2017, width=7, height=7)
directory_2019 <- ggmap(roadmap, base_layer = ggplot(subset(directory, year==2019), aes(lon, lat))) +
  coord_cartesian() +
  geom_point(alpha=0.1, stroke=0, position=position_jitter(0.02, 0.02)) +
  geom_density2d()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
xdensity_2019 <- ggplot(subset(directory, year==2019), aes(x=lon)) +
  geom_histogram() +
  scale_x_continuous(limits=c(-126, -113)) +
  labs(y=element_blank(), x=element_blank())

ydensity_2019 <- ggplot(subset(directory, year==2019), aes(x=lat)) +
  geom_histogram() + 
  coord_flip() +
  scale_x_continuous(limits=c(31.5, 41.5)) +
  labs(y=element_blank(), x=element_blank())

distribution_2019 <- gridExtra::arrangeGrob(xdensity_2019, blankPlot, directory_2019, ydensity_2019, 
                                            ncol=2, nrow=2, widths=c(4, 1), heights=c(1, 4))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 122 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 123 rows containing non-finite values (stat_density2d).
## Warning: Removed 123 rows containing missing values (geom_point).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 119 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
grid::grid.draw(distribution_2019)

ggsave("../images/distribution_2019.png", distribution_2019, width=7, height=7)

Final plots - usda progression

usda_2015 <- ggmap(roadmap, base_layer = ggplot(subset(usda, year==2015), aes(lon, lat))) +
  coord_cartesian() +
  geom_point(alpha=0.2, stroke=0.5, position=position_jitter(0.1, 0.1)) +
  geom_density2d()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
usda_xdensity_2015 <- ggplot(subset(usda, year==2015), aes(x=lon)) +
  geom_histogram() +
  scale_x_continuous(limits=c(-126, -113)) +
  labs(y=element_blank(), x=element_blank())

usda_ydensity_2015 <- ggplot(subset(usda, year==2015), aes(x=lat)) +
  geom_histogram() + 
  coord_flip() +
  scale_x_continuous(limits=c(31.5, 41.5)) +
  labs(y=element_blank(), x=element_blank())

usda_2015 <- gridExtra::arrangeGrob(usda_xdensity_2015, blankPlot, usda_2015, usda_ydensity_2015, 
                                            ncol=2, nrow=2, widths=c(4, 1), heights=c(1, 4))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 2 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 2 rows containing non-finite values (stat_density2d).
## Warning: Removed 2 rows containing missing values (geom_point).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 1 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
grid::grid.draw(usda_2015)

ggsave("../images/usda_2015.png", usda_2015, width=7, height=7)
usda_2017 <- ggmap(roadmap, base_layer = ggplot(subset(usda, year==2017), aes(lon, lat))) +
  coord_cartesian() +
  geom_point(alpha=0.2, stroke=0.5, position=position_jitter(0.1, 0.1)) +
  geom_density2d()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
usda_xdensity_2017 <- ggplot(subset(usda, year==2017), aes(x=lon)) +
  geom_histogram() +
  scale_x_continuous(limits=c(-126, -113)) +
  labs(y=element_blank(), x=element_blank())

usda_ydensity_2017 <- ggplot(subset(usda, year==2017), aes(x=lat)) +
  geom_histogram() + 
  coord_flip() +
  scale_x_continuous(limits=c(31.5, 41.5)) +
  labs(y=element_blank(), x=element_blank())

usda_2017 <- gridExtra::arrangeGrob(usda_xdensity_2017, blankPlot, usda_2017, usda_ydensity_2017, 
                                            ncol=2, nrow=2, widths=c(4, 1), heights=c(1, 4))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 2 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 2 rows containing non-finite values (stat_density2d).
## Warning: Removed 2 rows containing missing values (geom_point).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 1 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
grid::grid.draw(usda_2017)

ggsave("../images/usda_2017.png", usda_2017, width=7, height=7)
usda_2019 <- ggmap(roadmap, base_layer = ggplot(subset(usda, year==2019), aes(lon, lat))) +
  coord_cartesian() +
  geom_point(alpha=0.2, stroke=0.5, position=position_jitter(0.1, 0.1)) +
  geom_density2d()
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
usda_xdensity_2019 <- ggplot(subset(usda, year==2019), aes(x=lon)) +
  geom_histogram() +
  scale_x_continuous(limits=c(-126, -113)) +
  labs(y=element_blank(), x=element_blank())

usda_ydensity_2019 <- ggplot(subset(usda, year==2019), aes(x=lat)) +
  geom_histogram() + 
  coord_flip() +
  scale_x_continuous(limits=c(31.5, 41.5)) +
  labs(y=element_blank(), x=element_blank())

usda_2019 <- gridExtra::arrangeGrob(usda_xdensity_2019, blankPlot, usda_2019, usda_ydensity_2019, 
                                            ncol=2, nrow=2, widths=c(4, 1), heights=c(1, 4))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 3 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 3 rows containing non-finite values (stat_density2d).
## Warning: Removed 3 rows containing missing values (geom_point).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 2 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
grid::grid.draw(usda_2019)

ggsave("../images/usda_2019.png", usda_2019, width=7, height=7)

Final plots - Napa zoom

directory_napa <- ggmap(napa_roadmap, base_layer = ggplot(subset(directory, year==2019), aes(lon, lat))) +
  geom_point(alpha=0.3, stroke=0, position=position_jitter(0.05, 0.05)) +
  geom_density2d()
directory_napa
## Warning: Removed 1909 rows containing non-finite values (stat_density2d).
## Warning: Removed 1921 rows containing missing values (geom_point).

ggsave("../images/directory_napa.png", directory_napa, width=7, height=7)
## Warning: Removed 1909 rows containing non-finite values (stat_density2d).

## Warning: Removed 1921 rows containing missing values (geom_point).
usda_napa <- ggmap(napa_roadmap, base_layer = ggplot(subset(usda, year==2019), aes(lon, lat))) +
  geom_point(alpha=0.5, position=position_jitter(0.05, 0.05)) +
  geom_density2d()
usda_napa
## Warning: Removed 151 rows containing non-finite values (stat_density2d).
## Warning: Removed 153 rows containing missing values (geom_point).

ggsave("../images/usda_napa.png", usda_napa, width=7, height=7)
## Warning: Removed 151 rows containing non-finite values (stat_density2d).

## Warning: Removed 153 rows containing missing values (geom_point).